/*

  Simple cross-product, outputs in correct triangular form.

  - Coalescing memory access in the Column reads
  - No memory coalescing in the Row reads
  - 50% coalescing in the write (float8s)
  - All threads read in all data

 */

#define NSTATION 512
#define NTIME 100
#define NPOL 2

__kernel void SimpleTriangle(__global const float2 *array_d, __global float2 *product_d)
{

  // get global thread ID
  int Col = get_global_id(0);
  int Row = get_global_id(1);

  if (Col > Row) return;


  //instantiate sum variables
  float2 sumXX = (float2)(0.0f, 0.0f);
  float2 sumXY = (float2)(0.0f, 0.0f);
  float2 sumYX = (float2)(0.0f, 0.0f);
  float2 sumYY = (float2)(0.0f, 0.0f);

  for(int t=0; t<NTIME; t++){
    float4 blurbCol = ((float4*)array_d)[t*NSTATION + Col];
    float4 blurbRow = ((float4*)array_d)[t*NSTATION + Row];
    float2 rowX = (float2)(blurbRow.x, blurbRow.y);
    float2 rowY = (float2)(blurbRow.z, blurbRow.w);
    float2 colX = (float2)(blurbCol.x, blurbCol.y);
    float2 colY = (float2)(blurbCol.z, blurbCol.w);

    sumXX.x += rowX.x*colX.x + rowX.y*colX.y;
    sumXX.y += rowX.y*colX.x - rowX.x*colX.y;
    sumXY.x += rowX.x*colY.x + rowX.y*colY.y;
    sumXY.y += rowX.y*colY.x - rowX.x*colY.y;
    sumYX.x += rowY.x*colX.x + rowX.y*colX.y;
    sumYX.y += rowY.y*colX.x - rowX.x*colX.y;
    sumYY.x += rowY.x*colY.x + rowY.y*colY.y;
    sumYY.y += rowY.y*colY.x - rowY.x*colY.y;
  } 

  ((float4*)product_d)[(Row*(Row+1)/2)*NPOL + Col*NPOL + 0] = 
    (float4)(sumXX.x, sumXX.y, sumXY.x, sumXY.y);
  ((float4*)product_d)[(Row*(Row+1)/2)*NPOL + Col*NPOL + 1] = 
    (float4)(sumYX.x, sumYX.y, sumYY.x, sumYY.y);

}
